NASA Contractor Report 178111 
ICASE REPORT NO. 86-39 


n)Asa-C£- 1 1/ 


ICASE 


NASA-CR-1781H 

19860020625 


INCIPIENT TRANSITION PHENOMENA IN 
COMPRESSIBLE FLOWS OVER A FLAT PLATE 


Gordon Erlebacher 
M. Yousuff Hussaini 


Contract Nos. NAS1-17070, NAS1-18107 
June 1986 


INSTITUTE FOR COMPUTER APPLICATIONS IN SCIENCE AND ENGINEERING 
NASA Langley Research Center, Hampton, Virginia 23665 

Operated by the Universities Space Research Association 


NASA 

National Aeronautics and 
Space Administration 

Langley Research Center 

Hampton, Virginia 23665 


LIBRARY 


AuO 9 


2 6 



bdG 



LANGLEY RESEARCH CENTER 
LIBRARY, NASA 
HAMPTON, VIRGINIA 



Incipient Transition Phenomena in Compressible Flows over a Flat Plate 


G.Erlebacher 

NASA Langley Research Center. Hampton, VA 


M.Y. Hussaini 

Institute for Computer Applications in Science and Engineering 


ABSTRACT 


The full three-dimensional time-dependent compressible Navier-Stokes equations are solved by a Fourier- 
Chebyshev method to study the stability of compressible flows over a flat plate. After the code is validated in 
the linear regime, it is applied to study the existence of the secondary instability mechanism in the supersonic 
regime. 
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1. Introduction 


Although the predominant interest in the stability and transition to turbulence lies in compressible flows, 
research undertaken over the past 15 years has mainly focused on incompressible phenomena because of their 
comparative simplicity. As a result, a fairly extensive collection of theoretical [1,2], experimental [3,4] and 
numerical [5,6] data has been accumulated and cross-correlated for incompressible channel and Blasius 
boundary-layer flows. Detailed studies of several instability mechanisms has led to a fairly comprehensive pic- 
ture of the incipient stages of transition in incompressible flows. One such mechanism involves the interaction 
of a two-dimensional primary unstable wave with two skewed waves. This secondary instability leads to what is 
commonly referred to as K-type breakdown whose signature is a peak-valley vortical structure [1]. These struc- 
tures have been successfully simulated numerically for incompressible boundary-layer flows [5,6]. Moreover, 
results from these simulations agree very well with theoretical and experimental findings. 

As a step towards a better understanding of the instability mechanisms in compressible flows, the K-type 
breakdown of laminar flow at high Mach numbers is studied in this report To this end, a three-dimensional, 
fully-spectral compressible Navier-Stokes code capable of direct simulation of parallel boundary-layer flows over 
a flat plate has been developed. The code is first validated in the linear regime against the unstable eigenfunc- 
tions of the compressible linear stability eigenvalue problem. Then the temporal evolution of a triad of waves 
superimposed on a parallel boundary-layer is followed up to the incipient breakdown. 

2. Equations 

The full unsteady, three-dimensional, compressible Navier-Stokes equations expressed in non-dimensional 
form are 

■^+V(pv) = 0, (1) 

-^ + V- (pv*v) = -Vp + (2) 

-|g- + v>- Vp = -upV- v>+ - — V- (|iVT) + O + F p ( 3 ) 

and 


yMip = pRT 


(4) 
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where v* is the velocity vector (u,v,w), p the density, T the temperature, p the pressure, ¥ the viscous stress- 
tensor 

Y=n (Vv*+ Vv* 7 ) + XV • V* % (5) 

and O is the dissipation function 


O = (Vv^+ Vv> r ) 2 + X(V ■ v*) 2 . 


( 6 ) 


The first and second viscosity coefficients p and X are related by the assumption of zero bulk viscosity 

X + |n = 0. (7) 

= (Fp u » F pv > 0) and F p are source terms that are described later. Variables have been non-dimensionalized 
as follows: distance with respect to boundary-layer displacement thickness 8*, velocity, temperature and density 
with respect to their free-stream values 7L, pL and pressure with respect to the dynamic pressure pLLL 2 . 
Viscosity and conductivity are scaled with respect to their free-stream values Ju* and kL. With these definitions, 
the Reynolds number becomes 


„ P- 8* U„ 
R e = — 


(B) 


In this report the Prandtl number 



( 9 ) 


K 

is assumed constant and equal to 0.72. C p is the specific heat at constant pressure. We assume a temperature 
dependent viscosity given by the non-dimensional form of Sutherland's law 


vcn = 


d + 


196.8 


(T + 


196.8 


2 

r 2 . 


( 10 ) 


In incompressible boundary-layer flows, it has been observed that the evolution to transition of ribbon- 
excited Tollmien-Schlichting (TS) waves occupies a space of approximately ten wavelengths [3]. An adequate 
resolution of the boundary-layer in all three directions over a ten wavelength streamwise distance would exceed 
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the present capability of today’s most powerful computers. In incompressible simulations this difficulty has been 
circumvented by the assumption of spatial periodicity in the two directions parallel to the plate [5], allowing the 
set up of a properly posed initial-value problem. This assumption allows the computational domain to be res- 
tricted to one TS wavelength in both the streamwise and spanwise directions. A similar periodicity assumption 
is made in the present work to make the compresssible simulations computationally tractable. In order that the 
mean parallel boundary-layer flow be a stationary solution to the Navier-Stokes equations, source terms must be 
added to equations (l)-(3). The forcing terms are given by 


F = — - 

Fp “ R t dz ^ dz } 


F - — - 

Fpv R t dz ^ dz ) 


Fp = 


1 a , dT n ^ i J U m \ z dV m \ 2i 

TJM. Tz ^ + + ( — } 7 


( 11 ) 

( 12 ) 

(13) 


In obtaining these relations mean pressure is taken constant across the boundary-layer, the mean velocity has 
non-zero streamwise (U m (z)) and spanwise (V m (z)) components, and the mean temperature is 7 m (z) . In what fol- 
lows, V m is set to zero. 


3. Algorithm 

The system of equations (1-4) is solved in conservative form under the parallel flow assumption. Periodi- 
city boundary-conditions in the streamwise and spanwise directions permit a Fourier representation of the primi- 
tive variables u, v, w, p and p. For example, the double Fourier decomposition of u is 


N x N 
—-1 -Z-i 
2 2 

N x N 

T*2 


,i (majc+nj3>) 


(14) 


where N x and N y are the total number of nodes in x and y directions. The periods of the physical domain in the 


streamwise and spanwise directions ( L x and L y ) are related to the wave numbers a and P by Ly=-p * 


2tz 


The code has the option to evaluate normal spatial derivatives either by finite-difference approximation or by a 
Chebyshev collocation method, in which case u^z 9 t) has the series representation 
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N, 

“<>(*,0 = x U ijk (t) T k (z ) 

tel 

where T k (z) is the Chebyshev polynomial of order k. The number of nodes in the normal direction is N z . In 
the absence of discontinuities in the solution, spectral collocation methods are far more accurate than typical 
finite-difference methods for a specified distribution of a fixed number of nodes. The physical domain (0, z^) 
in the direction normal to the plate is mapped onto the computational domain with an algebraic mapping which 
takes the form 

_ *1/2 z max (1 + t]) 

*max (*max 2 Zj/j) 

for the spectral representation in the normal direction and 

z max *1/2 £ 

2 rr 1 

(*max “ *1 a) + % (2*1/2 ~ *tok) 

when finite-differences are used. In the latter case, the computational variables £, ranges from 0 to 1. In terms 
of the standard Chebyshev spectral node distribution is given by 

T| = — COS(TC^) , 

One half the normal nodes are located between z=0 and z=z m . 

In all the mean flow calculations considered thus far, the critical layer is always found to lie between 5* 
and 1.55* from the wall and special care is excercised to resolve it. Setting z V2 = 2 is found to be optimal for 
resolving the relevant flow structure in the chosen parameter range. 

Currently, the code is fully explicit. A third order low-storage Runge-Kutta method [8] is used for time 
discretization. For the equation 

U t = F, 

H X =F 0 

U\ = U 0 + -j Hi 


such a scheme leads to 
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h 2 = &f 1 -1h 1 

U 2=Ui + ^H 2 

H 3 = 8//?2 "l§ //2 
u i = U 2 + ±H y 

which only requires two storage arrays per variable. Theoretically the scheme is unconditionally stable for a 
CFL below 0.55; however it is empirically found that above a CFL of 0.2 that the algorithm is unstable. The 
time step is thus limited by 

At < 0.2 min [— + ^ + n — ] . 

grid U V W 


4. Initial Conditions 

The initial conditions consist of a triad of waves superimposed on a mean flow. Mean flow profiles are 
generated from the solution to the similar compressible boundary-layer equations with zero pressure gradient 
and zero heat transfer at the wall. In non-dimensional form, these are [9] 

#" + (Hp/")' = 0, 

— ( p p T') '+fT' + (7-1) Mi p p T,/" 2 = 0 , 

pT=l, 

with boundary conditions 


m =/'( 0 ) = 0 , 
/'(-)= 1, 
7X°°)=1, 

T '(0) = 0 . 

A prime denotes differentiation with respect to the similarity variable 


H = c(M w TJ) 


_Z 

>/2 


U~ &*(x) p„ 
M- 


Z_ 

& 


R t (x) 


(15) 
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where Z is related to the normal coordinate z by the transformation 

Z 

Z = | p dz 

The coefficient c in (15) is a function of free-stream temperature and Mach number. 

The stability of flows is studied at a prescribed supercritical Reynolds number R t (x 0 ). After freezing the 
streamwise dependence of the mean flow profiles at x=Xq, they are extended ova* the entire plate. Thus, the dis- 
placement thickness is constant, and the mean flow is a function only of the coordinate normal to the plate. The 
normal velocity component is set to zero, and U m =f (T|). 

In the simulation which is to be described, the mean flow is initially perturbed by a 2-D Tollmien- 
Schlichting wave and two 3-D waves. For example, when perturbed by a single wave, the initial streamwise 
velocity is the real part of 

u = t/ m + i5(z)e‘ (cu + p >" <0 ' ) 

for real streamwise and spanwise wavenumbers a and p. Similar relations hold for the remaining variables. The 
complex growth rate co is calculated by an eigenvalue code [10]. From the real part of co, the period of the wave 
can be determined. Growth or decay of the wave depends on the sign of co, , the imaginary part of co. A stable 
wave corresponds to co 4 - negative, while a positive co, indicates a growing mode. 

5. Boundary conditions 

Under the assumption of parallel flow, all the variables are periodic in the streamwise and spanwise direc- 
tions. No slip conditions are applied to the velocities at the wall which is adiabatic. In the far-field, all the 
variables are frozen at their initial values. 

6. Results 

The coordinate system and nomenclature used for the flat plate geometry is illustrated in figure 1. Com- 
plex eigenfunctions and eigenvectors are obtained by numerically solving the eigenvalue problem associated 
with the linear stability of compressible parallel flows [10]. Table I summarizes the input parameters. Mean 
streamwise velocity and mean temperature profiles at Mach 0.5 and 4.5 are compared with each other in figures 
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2-3. As seen from figure 2, at Mach 4.5 the mean temperature across the boundary-layer varies by more than a 
factor of 4. This is in sharp contrast with the flat temperature profiles at Mach 0.5 (fig. 3). The complex eigen- 
functions are normalized to insure a maximum streamwise perturbation velocity of unity. Streamwise velocity 
profiles are illustrated in figures 4-5. With the normalization imposed on the eigenfunctions, it is immediately 
apparent from figure 6 that the extrema of the temperature variation are much more severe at the higher Mach 
number. 

Added to the mean flow, the computed three-dimensional eigenfunctions serve as initial flow-field data 
for the full Navier-Stokes code. For example, the initial streamwise velocity distribution is 

i/(z) = U m + e 2D cos(0(z)+ocx+Py) 

where 0(z) is the phase of the streamwise velocity eigenfunction. Similar relations are satisfied by the remaining 
velocity components, pressure and density. A value of .001 for e 2 d is sufficient to insure the absence of non- 
linear interactions of the fundamental mode with its higher harmonics. Figures 7-8 summarize the results of the 
linear test. Plotted in these figures is log{E K ) versus time, where E K is the perturbed kinetic energy 

Ek = J j Pm [(u-UJ 2 + v 2 + w 2 ] dx dy dz . 

Comparisons of the growth rates predicted by the Navier-Stokes code against linear results are made using both 
Chebyshev collocation and finite-differences in the direction normal to the wall. At Mach 0.5, the growth 
curves of perturbed kinetic energy for a 33 collocation point normal distribution and a distribution with 129 
finite-difference nodes are almost indistinguishable from each other and from the growth predicted by the linear 
eigenvalue code. However, 65 Chebyshev collocation nodes are required to match the predicted linear growth 
curves at Mach 4.5. In this case, 129 finite-difference nodes are not sufficient to resolve the structure of the 
eigenfunctions near the wall. The need for extra resolution at the higher Mach number is explained by the more 
complicated structure of the eigenfunctions presented in figures 5 and 6b. At Mach 4.5, the displacement thick- 
ness is almost an order of magnitude greater than at Mach 0.5, but in units of 5*, the distance of the critical 
layer from the wall always lies between 1 and 1.5. 

For the non-linear simulation, the parameters are: AC=4.5, *=10000, 7L=110.85°i?, a= 6 and p=1.03923. 

The spanwise wavenumber is chosen to maximize the growth rate of the 3-D mode in order to accelerate the 
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onset of the instability. The initial wave angle of the three-dimensional wave is 60° . 

The amplitudes of the 2-D and 3-D waves are respectively £20 and £jd ■ Initially, the flow variables are 

set to 


u(x,y,z) = U m + 1 2 D (z) cos (8^(2) + ou + Py) + 

twUiD /cos(0„ JO (z) + ax + Py) + cos(0„ w (z) + ax - Py)7 

v(x,y,z) = 620(2) V20 cos(0 V2D (z) + ax + Py) + 

£3 d v 3 d [cos (0 VjD (z) + ax + Py) - cos(0 Vjo (z) + ax - Py)/ 

w(x,y,z) = £20 m> 2 d{z) cos(0h Ijd (z) + ax' + Py) + 

63 d W 3D [ COS (0 Wjo (z) + ax + Py) + cos(0 WjD (z) + ax - Py)/ 

p(x,y,z) = p m +£ w Pi D (z) cos(0 P2O (z) + ax + Py) + 

e 3 o p3Z) l cos(0 PjD (z) + ax + Py) + cos(9 Pjd (z) + ax - py)7 


p(x,y,z) =p m + £2D P2d(z) cos(Q Pid (z) + ax + Py) + 

£30 PiD [ cos (0^ o (z) + ax + Py) + cos(0 Pjo (z) + ax - 3y)7 


Starting amplitudes are z 2 &=-054 and £j/>=.012. The 2-D growth rate predicted by linear theory is 
©2/r=.5011+.00203i with a period of T=12.54. The 3-D growth rate, coj^.09765+,01098/, is substantially larger 
than its 2-D counterpart, a property typical of compressible flows as the Mach number increases. Above Mach 
3, the phase angle that produces maximum growth is between 55° and 65° [2]. As stated earlier, the primitive 
variables are expanded in Fourier series in the x and y directions. The (i,j) mode of u is defined by the 
coefficient Uq in (14). Figure 9 shows the evolution of the kinetic energy of selected modes as a function of 
time. The energy content of mode (ij) is defined by 


Eij ~ logio ( J- 


Ujj- + v,/ + V 


dz). 


The energy content of the (1,0) mode remains fairly constant The spanwise modes (0,1) and (0,2) continuously 
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grow until their strength is apparently sufficient to trigger the growth of the (1,1) mode, which is largely respon- 
sible for the presence or absence of the secondary instability. Prior to the rapid growth of the (1,1) mode, it 
goes through a stable interval starting at about 3 periods. Once the (0,1) mode is within 0.5 decades of the fun- 
damental 2-D mode, the (1,1) mode suddenly becomes unstable again. At this stage, non-linear interactions 
between the 2-D and the 3-D mode begin to dominate the evolution of the flow. 

Spanwise vorticity contours are illustrated in figure 10 in successive spanwise planes (normal to the span- 
wise direction) in the peak and adjacent planes. The initial formation of the high-shear layer after 5.5 periods 
(computed on a 36x16x65 grid) is apparent in this figure. The downward curving of the shear layer near its 
maximum is similar to the structure found in incompressible flow in the early stages of transition [5,6]. Stream- 
wise and spanwise vorticity contours in streamwise planes (the normal to the plane is in the streamwise direc- 
tion) are shown in figures 11-12 after 5.5 periods in planes at 0.0, 0.25, 0.5 and 0.75 wavelengths along the 
streamwise direction. Saddle points in the spanwise vorticity correspond to maximum streamwise vorticity. If 
these maxima are followed in figures 12a-d, the uplifting of the vortex tube becomes apparent, although it is 
rather flat 

7. Conclusions 

A fully spectral, 3-D, compressible Navier-Stokes code specialized to parallel flow over a flat plate has 
been presented. It has been demonstrated that in the linear regime an initial eigenfunction grows at the same rate 
as predicted by a linear eigenvalue code in both subsonic and supersonic flow. As a test of the non-linear 
behaviour of the code at high Mach numbers, initial conditions are set up to numerically generate, for the first 
time, a high Mach number peak-valley vortical structure similar to that observed and computed for incompressi- 
ble flows. 
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9. Table 


Case 

Mach nb. 

a 

P 

R< 

(°R) 

I 

0.5 

.24933 

.20944 

1200 

520.00 

II 

4.5 

.4652 

.8057 

10000 

110.85 


Table I: parameters input to the linear eigenvalue code 
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11. Figure Captions 

Figure 1: Parallel mean flow over a flat plate. The coordinate axes are labelled according to accepted 

nomenclature. 

Figure 2: Mean profiles U m (solid line) and T m (dashed line) at Mach 0.5 (see table I) along the direction 

normal to the plate. 

Figure 3: Mean profiles U m (solid line) and T m (dashed line) at Mach 0.5 (see table I) along the direction 

normal to the plate. 

Figure 4: Mach 0.5 complex eigenfunctions of the linearized compressible Navier-Stokes equations at 

i?*=1100 along the direction normal to the plate: u 2D (4a), u 3D (4b), T 2D (4c) and T 3D (4d). 

Figure 5: Mach 4.5 complex eigenfunctions of the linearized compressible Navier-Stokes equations at 

R e =10000 along the direction normal to the plate: u 2D (5a), u 3D (5b), T 2D (5c) and T 3D (5d). 

Figure 6: Comparison of temperature eigenfunction profiles at Mach 0.5 (6a) and Mach 4.5 (6b). 

Figure 7: Growth curves of perturbed kinetic energy at Mach 0.5 (case I) for 1/2 period with 33 collocation 

nodes (x), 33 (+), 65 (*) and 129 (x) finite-difference nodes. The 129 finite-difference, 33 collo- 
cation and the theoretical curve are indistinguishable from each other. 

Figure 8: Growth curves of perturbed kinetic energy at Mach 4.5 (case II) for 1 period with 129 finite- 

difference nodes (*), and 33 (x) and 65 (O) collocation nodes. The theoretical curve is indistin- 
guishable from the 65 spectral node curve. 

Figure 9: Evolution of the kinetic energy of selected Fourier modes. The energy is plotted on a logarithmic 

scale. 

Figure 10: Spanwise vorticity in successive spanwise planes after 5.5 periods (^=4.5,^=10000.0). The 

planes are located at x=0 (10a), LJ4 (10b), and LJ 2 (10c) where L x is the streamwise period. 
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Figure 11: 


Figure 12: 


Figure 10c is the peak plane. The planar cuts extend to 2=2 in the normal direction. 

Spanwise vorticity in successive streamwise planes after 5.5 TS periods (^=4.5,^=11000.0), 
The planes are located at y=0 (11a), L/4 (lib), Zy 2 (11c), and 3Zy4 (lid) where L y is the span- 
wise period. The planar cuts extend to z=2 in the normal direction. 

Streamwise vorticity in successive streamwise planes after 5.5 TS periods (^=4.5^=10000.0). 
The planes are located at y~0 (12a), Zy4 (12b), Zy 2 (12c), and 3Zy4 (12d) where L y is the span- 
wise period. The planar cuts extend to z=2 in the normal direction, 


\ 



z (normal) 


14 



Figure 















KINETIC ENERGY 


17 




Figure 7 Figure 8 


— 3 . 00 
- 3 . 40 
— 3 . 80 
— 4 . 20 
- 4 . 60 
— 5 . 00 
- 5 . 40 
- 5 . 80 
— 6 . 20 
- 6 . 60 
- 7 . 00 

0 . 00 1 . 20 2 . 40 3 . 60 4 . 80 6 . 00 

PERIODS 



Figure 9 















Standard Bibliographic Page 


1. Report No. NASA CR- 178 111 
ICASE Report No. 86-39 


2. Government Accession No. 


3. Recipient’s Catalog No. 


4. Title and Subtitle 

INCIPIENT TRANSITION PHENOMENA IN 
COMPRESSIBLE FLOWS OVER A FLAT PLATE 


5. Report Date 

June 1986 


6. Performing Organization Code 


7. Author(s) 

Gordon Erlebacher 


and M. Yousuff Hussaini 


8. Performing Organization Report No. 

86-39 


9. Performing Organization Name and Address 

Institute for Computer Applications in Science 
and Engineering 

Mail Stop 132C, NASA Langley Research Center 
Hampton. VA 23665-5225 


10. Work Unit No. 


11. Contract or Grant No. 

NAS1-17070, NAS1-18107 


12. Sponsoring Agency Name and Address 

National Aeronautics and Space Administration 
Washington, D.C. 20546 


13. Type of Report and Period Covered 
Report- 


14. Sponsoring Agency Code 

505 31 83 01 


15. Supplementary Notes 

Langley Technical Monitor: 
J. C. South 

Final Report 


Submitted to Proc. 10th Inter. 
Conf. on Num. Methods in Fluid 
Dynamics 


16. Abstract 

The full three-dimensional time-dependent compressible Navier-Stokes 
equations are solved by a Fourier-Chebychev method to study the stability of 
compressible flows over a flat plate. After the code is validated in the 
linear regime, it is applied to study the existence of the secondary 
instability mechanism in the supersonic regime. 


17. Key Words (Suggested by Authors(s)) 

compressible Navier-Stokes equations, 
boundary layer, secondary instability 


18. Distribution Statement 

34 - Fluid Mechanics and Heat 
Transfer 

Unclassified - unlimited 


19. Security Classif.(of this report) 

20. Security CIassif,(of this page) 

21. No. of Pages 

22. Price 

Unclassified 

Unclassified 

21 

A02 


For sale by the National Technical Information Service, Springfield, Virginia 22161 

NASA-Langley, 1986 










End of Document 



